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SUWARY 

A review of past parabolized Navler-Stokes applications Is presented. The equations, boundary 
conditions, the numerical method and the grid generation are all discussed. Results ranging from the low 
supersonic regime to the hypersonic regime are Included. 


1. INTRODUCTION 

The Intent of this long-term project has been to develop a CFD tool that will be capable of analyzing 
the viscous supersonic/hypersonic flow about realistic configurations. To meet this challenge, a tech- 
nique was developed that Is capable of predicting the flow in regions of canopies, wings, and canards In 
addition to the usual simple symmetric configurations. The technique also had to allow for Interactions 
between aerodynamic surfaces such as the vortex interaction between canards and wings. 

The NASA Ames Parabolized Navler-Stokes (PNS) code was used as the mainline procedure to simulate 
numerically the viscous supersonic flow over these generic configurations. The parabolized approximation 
to the Navler-Stokes equations assumes that the flow Is supersonic In the streamwlse direction, and that 
the subsonic flow In the viscous sublayer Is always positive in the streamwlse direction. Thus, flows 
with large streamwlse separation and flow reversals are excluded from treatment under the foregoing 
assumptions. However, crossflow separations are permitted. 

Under these assumptions, the Navler-Stokes equations become parabolic In the streamwlse direction, 
enabling a marching solution procedure which Is computationally desirable and efficient. The form that Is 
presented here was developed by Schlff and Steger (Ref. 1). It differs from the first f Inlte-dlfference 
Implicit marching algorithm in delta form of Vlgneron et al. (Ref. 2) In the way in which the streamwlse 
convection flux vector Is treated. The turbulence model that has been used in this project Is the 
Baldwln-Lomax model (Ref. 3). In some cases the turbulence model was modified to allow for a better 
comparison with experiment In the supersonic range. The boundary conditions are the usual viscous no-slip 
at the wall and a characteristic procedure is used to fit the bow-shock wave which Is the outermost bound- 
ary. Since the equations are cast In conservation-law form, all discontinuities within the flow domain 
are predicted accurately. Owing to the possible complex nature of the configuration at each axial loca- 
tion as one proceeds down the body, an elliptic grid generator Is employed to discretize the flow domain. 
In addition, an equilibrium air capability has been incorporated into the code. It uses the curve fits of 
Tannehill et al. (Ref. 4). Body shapes that have been marched over by the present algorithm Include 
pointed cones (Ref. 5), sphere-cones (Ref. 6), three-dimensional reentry vehicles (Refs. 6-8), an oglve- 
cylinder-boattall (Ref. 9), a delta-wing configuration (Refs. 2, 10), the X-24 lifting body (Ref. 7), the 
Space Shuttle (Ref. 11), finned configurations (Refs. 5, 12), and, most recently, a generic fighter con- 
figuration (Refs. 13, 14). A representative number of these results are included in the present work. 


2. GOVERNING EQUATIONS 

The governing equations in the two base coordinate systems (Fig. I) are transformed under a general- 
ized coordinate mapping 

t = C(a) 

n = n(a,b,c) (1) 

c = c(a,b,c) 

where the system (a,b,c) is either the Cartesian or cylindrical coordinate. The surface c * 0 is the 
body surface, and a c = constant surface is the outer boundary or the bow-shock surface. Since the 
solution is marched along the axial direction, a, the c = constant surface is chosen to be the axis- 
normal crossflow plane; therefore, c is not a function of b or c. 

For the Cartesian form a, b, and c are more commonly known as x, y, and z. In many cases, the 
cylindrical form provides a better choice of governing equations than the Cartesian for the particular 
geometry of concern. And, for that reason, this report and the code contain both the Cartesian and cylin- 
drical form of the Navler-Stokes equations. In the cylindrical form a, b, and c are cylindrical coordi- 
nates x, e, and r. Where r is the length of a ray extending from the x-axis in a constant x plane, 
and 8 Is the angle of that ray with the y-axis (z = 0) corresponding to the constant x plane. 

In the thin-layer approximation, viscous diffusion in the axial direction as well as along the sur- 
face of the body in the crossflow plane (n-direction) is neglected. The thin-layer approximation used 


1 



here as a starting point is not essential to the development of the PNS equations. Tests of high- 
Reynolds-number flow In two dimensions (Re = 10 4 ) have always shown that the additional terms left out 
because of the thin-layer approximation cannot be detected to four-place accuracy in a typical finite- 
difference solution. Finally, the approximation termed the "parabolized approximation" is used. The 
diffusion terms in the marching direction, t, are neglected. A relatively stable marching procedure 
results if the following conditions are satisfied: 1) the outer inviscid flow region is supersonic rela- 
tive to the downstream marching coordinate; 2) the viscous layer does not separate in the marching direc- 
tion; and 3) the sublayer model of Schiff and Steger (Ref. 1) is used. 


By invoking these approximations, the normalized Reynolds averaged Navler-Stokes equations transform 
under the generalized mapping as follows: 


aE t aF ' aG : _ 1 3G v 

at an a? s ' Re a; 



( 2 ) 


The vectors E, F, and G represent the inviscid flux vectors in the three transformed directions. 
They are written as, 
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where q is the vector of dependent variables. In these expressions, p and p represent the density and 
pressure, respectively. The variables v and w are the velocity components in the y and z directions 
for Cartesian coordinates, and In the e and r directions for cylindrical coordinates, u is the veloc- 
ity component in the x direction. The vectors S and S y in Eq. (2) are source terms because of the 
cylindrical coordinate formulation. When using Cartesian coordinates, these vectors do not exist. 

The Jacobian of the transformation, J, is related to the metrics by 


J = 5 a ( Vc ' Vb> + <b ( Va ' Vc* 

+ S C <Va - Vb> 

(4) 

where the metric quantities are given by 



'a = J < b n c c * W * °a = J <V>c ‘ Vc* 1 

‘a * J Vn * c S b n ) 


<b = J < c n a c - 4 n C c } * "b ‘ J <Vt - W • 

c b = °<Vn - a i C n> 

(5) 

^C = J(a n b c ' b nV * n C = J(b <; a c ; Vcl* 
The total energy, e, is defined as, 

, 'c. s J(, E b , ■ W 
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(6) 

where c is the internal energy. The contravariant velocities 

U, V, and W. are given as 


U = t a u + C b v ♦ C c w 


(7) 

V * HjU + n^v + n c w 


(8) 

W = c a u + V + < c w 


(9) 
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The only remaining viscous flux vector after the parabolized assumptions are taken can be written as, 
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where the quantities raj, m 2 , and m 3 are given as, 
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for these equations v, k, and T represent viscosity, thermal conductivity, and temperature. The second 
viscosity coefficient has been taken as x = -2/3u. In the Cartesian coordinate formulation, 

a = x 


b * y (14a) 

c = z 


and 



(14b) 


and In the cylindrical coordinate formulation, 

a = x 
b = 9 
c = r 


(14c) 


with 


and 


with 


where 


Is replaced with 
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Having defined the Cartesian and cylindrical coordinate formulations and the necessary relations In 
order to switch from one to the other, subsequent discussion will be carried out In the Cartesian coordi- 
nate formulation only. • 
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Nondimensional ization 


The nondimensional form of all the variables has been used throughout. They can be related to the 
dimensional form (denoted by the superscript -) by the following definitions. 


i p [ 

m co 


= a/a^ 

or co 

Re = — — 


a 2 ; 

= u/a^ 

Pr = 

T_E_ 

= 5/a„ 

T = T/T„ 
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v = u/u„ 
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= 

- —2 
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The subscript - refers to values of the freestream. The value of the dimensional length t is 
generally taken as one unit of length (unity), and thereby allows one to work with dimensional lengths 
which can present less of a problem when coding the geometry, or when reading output. The nondimensional 
quantity, Pr, Is a form of the Prandtl number and for a perfect gas is related to the Prandtl number by 
Pr = ( Y - 1 ) Pr . 


Closure 

The system of equations is closed with an appropriate equation of state and with the definition of 
properties. Real gas effects must not be neglected at high Mach numbers or for high heating cases. It Is 
In these cases that the perfect gas assumption and constant property assumptions are Invalid. Therefore, 
an appropriate equation of state must be used. A generalized equation of state can be written in the 
form, 

P = P(e.o) (17) 

or In a more specific form, 

P = (y - 1)pe (18) 

where 

y = y(e,p) (19) 

Here y Is the ratio of specific enthalpy to specific internal energy. It is only for a perfect gas 
that this ratio Is a constant and equal to the ratio of specific heats, y. 

The expression of the pressure as a function of Internal energy and density simplifies the calcula- 
tion because both internal energy and density can be directly expressed In terms of the dependent varia- 
bles. In contrast, an equation of state that is dependent on temperature requires more computation since 
temperature Is not a dependent variable and, in general, is not a readily known function of the dependent 
variables. 

The three relations needed to close the set of equations are the relation of temperature, viscosity, 
and thermal conductivity to the dependent variables. 


T = T(c,S) 

(20) 

V - v(e.o) 

(21) 

E = E(c,c) 

(22) 


Again the relation to e and p simplifies the calculation. 

A variable that Is needed for the calculation of Mach number is the speed of sound, a. The speed of 
sound Is a direct result of the governing equations, including the equation of state. For consistency, 
the speed of sound should be derived from the governing equations, and not approximated with a curve fit 
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or a table look-up of experimental data. An eigenvalue analysis reveals the speed of sound. r o: the 
generalized equation of state in £q. (17), the speed of sound is found to be 
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( 23 ) 


If the pressure derivatives in the calculation of the speed of sound are handled In the same manner 
as In the rest of the code, then consistency is assured. With this reasoning, the equation of state 
should not only return accurate values of pressure, but should also return accurate values of the first 
derivative of pressure. However, If the fluid is operating under the perfect gas equation of state, then 
the perfect gas relations 



are used. 


(24) 


3. BOUNDARY CONDITIONS 

Boundary conditions must be applied at: (1) an initial plane of data at c = constant, (2) both 
extremes of the computational domain in the c -coordinate, and (3) at two planes of data in the circumfer- 
ential (n) direction (unless periodic). 

An accurate initial data solution is needed. The parabolized Navier-Stokes prediction procedure can 
generate Its own appropriate starting solution for conical bodies, or it can be started from a time depen- 
dent blunt-body program. Both options are discussed subsequently. 

The boundary conditions specified at the surface of the body always include the no-slip and 
no-blowing conditions. Options are available that allow specifications of either an assigned wall tem- 
perature or an adiabatic wall condition. 

For symmetric bodies at zero angle of yaw, the synmetry conditions in the n = 0 and n = n max 
planes are Imposed. Otherwise, periodic boundary conditions are imposed in n . Although the periodicity 
condition requires more computer storage, it is the easiest boundary condition to Impose implicitly in 
the n-dlrection without stability problems. 

The outer boundary of the computational domain Is either specified as a known free stream or as the 
peripheral bow shock which Is determined as the solution is marched downstream. Shock-fitting logic is 
used In the latter case. 

The various boundary conditions are discussed in the following subsections: the Initial conditions, 

the body-boundary conditions, shock-fitting procedures, and the circumferential boundary conditions. 

Initial Data 

Consider a flowfield which starts at an initial plane defined by c? = constant and extends down- 
stream to larger values of t. Two planes of data are required to start the calculations, one at 
tj = constant and the second at ^ = constant (where tj < t ? ). In general, each of the initial planes 
will be divided into a grid defined by combinations of c,n coordinates'. Initial flowfield data given In 
the format a, pu, pv, pw, and e are to be supplied at each of the t.c.n grid points. Here, the 
Cartesian-velocity components are employed and the variables are not normalized by Jacoblans. It is not 
necessary that the Initial data be supplied on a grid which is exactly the same as the one to be used 
downstream, since the PNS code will automatically interpolate the Initial data onto an appropriate grid. 

There are always potential problems when numerical solutions are matched along a common boundary. 
Inconsistencies can be generated by differences in the governing equations, the solution procedures, the 
computational grids, or any number of other sources. A corononly encountered systematic error Is Intro- 
duced when the blunt-body procedure does not Include a constant pressure, viscous-sublayer model. Even 
when the prediction procedures are well matched, significant errors can be Introduced by the Interpolation 
from one grid to the other. With a little forethought and judgment, the careful investigator should be 
able to keep the effects of such inconsistencies well within the truncation error which Is always present 
In the downstream PNS solution. 

For body configurations at moderate angles of attack and yaw, and almost conical in shape near the 
nose, one can use the conical approximation to obtain starting data. Indeed, by employing the so-called 
step-back principle, one can use the PNS code to self-start the calculations. This principle Is based on 
the supersonic cone rule at zero angle of attack, which Is approximate for turbulent flow and exact for 
laminar flow. The essential idea is that the solution over the conical part of the body near the nose can 
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be approximated by the solution over an equivalent pointed cone. Therefore, the similarity principle 
appl ies, and' the solution Is' only' a function of' the' similarity coordinate’ and varies over the^downstream 
distance by a scaling factor’ only’; ’ The flow,’ be'lng'self-slmilar,- can therefore be : calculated’ by'ahVitera- 
tive technique. Using such an approach, one advances the data from t station j to j + 1 and then 
resets the flowfield solution back to j at an- n and ? station along the same conical ray. The flow 
variables are continuously readvanced to j + 1 and;,then'set back until a steady-state solution is 
reached. It Is understood that this Is only correct In the invlscld case. The slight error In the bound- 


ary layer thickness caused by this approximation In the viscous case is willingly accepted because of the 
computatTbnaieffi’ci^ncy 'involved; ■' ' ■') ." ! r - ■’ ■' ’ - a”**-* •< 


15 ■ J , Ax 1 symmetric; and’ three-dimensional 1 viscous, blunt-body solutions from.whlch to start the PNS program 
can be obtained using the- time-dependent, ’ blunt-body code developed by- NASA Ames . (Refs. 15, 16). i .The 
blunt-body program requires that the computational volume completely encompass; the embedded; subsonic 
region at the nose of the vehicles. At small angles of attack, i.e., o < 30 for many blunted vehicles, 
the entire embedded subsonic region will be on the spherically tipped portion of the nose. Such solutions 
can be rotated and used at other angles of attack, .provided that the embedded region remains on only the 
spherically tipped region of the nose. 


Body Boundary Condition 

For a viscous calculation, five boundary conditions must be specified at the body. The no-slip con- 
dition supplies three of these since all three velocity components must be zero at the wal 1. .^ Boundary- 
layer theory provides the fourth as the pressure derivative is zero in the normal direction to the body 
surface. The fifth boundary condition is the specification of either a wall temperature, or the specifi- 
cation of a heat flux. For the specified wall-temperature case, the boundary .conditions, can be^ written in 
the form 
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where the;, subscripts are associated .with , il . .and denote. the grid, point numbers in the -direction. 

■^iiall - 1s the. specif led wall temperature. Here a is,-0;or 1/3 for first-, .or,. second-order .accuracy.. The 
specif ied, heat if.lux case can. be expressed, as.- a,, specified .temperature, gradient at the. wall if .radiation 
effects are neglected. The velocity, pressure and temperature-grad1ent..boundary i cond1t1onS! can be.- 
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Again, the boundary condition Is corrected explicitly after an Implicit marching step Is completed. As 
long as the outer boundary Is a sufficient distance from the bow shock, this boundary condition is stable. 

The second outer boundary scheme Is shock fitting. In this method,. the outer boundary Is the bow 
shock. For this method, not only the fluid dynamic variables at the shock are calculated, but the varia- 
bles associated with shock propagation are calculated. A new vector of dependent variables Is defined as 

q = J"*(p,8U,pv,pw,e,y t ,z t ) (29) 

where y t and are the two unknown metric quantities that define the direction of shock propagation. 

For the seven unknowns at the shock, there must be seven equations. Five of the seven equations are 
the Ranklne-Hugonlot jump conditions across the shock which can be expressed In the form 

lG) j+1 = J' 1 (ljj +1 (E j+1 - EJ + c 3+1 (F J+1 - FJ + c 3+1 (G 3+1 - GJ) = 0 (30) 

The brackets denote a differencing across the shock, and the metric quantities are those at the point 
(J + l,k,l max ). 

In Eq. (30), not only are the flux vectors unknown at j + 1 but the metric quantities are also 
unknown, since they can be written In terms of the dependent variables y^ and z^. 

The fitted shock location at J + 1 Is found In three steps. The first step Is to predict a grid 
at j + 1. The shock slope at j Is extrapolated to define a predicted shock location at j + 1. A grid 
system Is then generated to span the gap between the body and the predicted shock. The second step is 
when the Implicit marching step from J to J + 1 is complete. The shock grid points at J + 1 are 
then freed and allowed to be Implicitly moved. A restriction on this movement results In the sixth equa- 
tion needed at the shock. 

The seventh and final equation needed to close the set of equations at the shock Is the compatibility 
equation. This equation is derived by first multiplying Eq. (2) by the Inverse of the Jacobian of the 
streamwlse convective flux vector and then by the left eigenvector corresponding to the only positive 
eigenvalue In the c-dlrectlon at the shock. This equation represents the flow field Influence on the 
Implicit shock solution. 

The compatibility equation combined with the grid point restriction and the Ranklne Hugonlot equa- 
tions closes the set of equations at the shock. The result Is a delta form of the shock equations with a 
7x7 matrix at the shock. 

Before a block tridiagonal Inversion can be done, the 7x7 matrix, D must be reduced to a 5x5 matrix. 
The shock equation can now be Inverted with the Interior flow-field equations to yield the Implicit result 
of the dependent variables at J + 1. 

An explicit correction of the shock variables. Is done to eliminate any linearization error that might 
have occurred In the Implicit marching procedure. In the real gas case, this Is done by accepting the 
Implicit result of density at the shock and iterating on an appropriate set of equations to yield the 
remaining unknowns at the shock. The unknowns are reduced to six since density Is known, and the compati- 
bility equation must be eliminated so that there are also six. The six remaining equations are the five 
Ranklne Hugonlot relations and the grid-point restriction equation. A Newton Iteration procedure Is done 
on these six equations until convergence of the shock dependent variables Is reached. 

Circumferential Boundary Conditions 

Boundary conditions must be Imposed at the bounding planes (k = 1 and k = k^) 1n the c1rcumfer en- 
tlal direction. Two types are available— periodic and symmetric (or reflective). Periodic-boundary con- 
ditions must be employed when a nonzero-yaw angle is specified. Symmetry conditions can be employed for 
zero-yaw angles for configurations that are symietrtc about a pitch plane. The use of a symmetry boundary 
will always result In much faster calculations owing to a fewer number of grid points, so this boundary 
condition Is used whenever It Is valid. 
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Periodic-Boundary Condition 

The periodic-boundary condition assumes that, in the circumferential direction, the first plane of 
data is adjacent to the last, i.e., the grid is assumed to wrap around the entire body. Relations of the 
form 

6 n q| k=l = \ (R - q k > (31) 

n 1 K max 

are employed to transmit this information to the solution procedure. 

Symmetry Plane Boundary Condition 

A second-order accurate, fully implicit set of reflective boundary conditions has been formulated and 
is included in the prediction procedure. 

This boundary condition assumes that the flow-field is symmetric (or reflective) about the pitch 
plane (y = 0). Thus, the flow-field properties at a given (x,y,z) point are related to the flow-field 
properties at the corresponding reflective point (x,-y,z) by 

y = -y (or n = -n) 

p(x.y.z) = p(x,-y,z) 

pu(x,y,z) = pu(x,-y,z) 

(32) 

pv(x.y.z) = -pv(x,-y,z) 

pw(x,y,z) = pw(x,-y,z) 

e(x,y,z) = e(x,-y,z) 

where the coordinate y and velocity component v change sign as they pass through the pitch plane. 

Currently, all symmetry planes occur for n = constant. Boundary-condition formulations are provided 
for the n derivatives of the F variables, one smoothing term and the turbulence model. 

4. NUMERICAL METHOD 

The algorithm used in the present work is based on the Beam-Warming delta form as applied by Schiff 
and Steger (Ref. 1) and other researchers (Ref. 2). It is first- or second-order accurate In the marching 
direction and second-order accurate in the spatial directions. A fourth-order dissipation term is 
appended to the algorithm as is an implicit smoothing term to further enhance the stability. 

The Implicit matching algorithm written in delta form without implicit smoothing is 

[A^ + (1 - a)A t (« n B j + 5| .C j - Re' 1 4 l « j )](q J+1 - q j ) 

- -(A| - A J’ 1 )^ + a (E 3 - E^- 1 ) - (1 - .)&' 

x <« !n 3 + 1 (E/J) j + ni + 1 (F/J) j + n 3 + 1 (G/J) j | + « U;j + 1 (E/J) j + c 3 + 1 (F/J) j + ^ + 1 (G/J) j ] - Re’ 1 ? r S J > 
n a y i. t a j . i. c 

- I(t x /J) j+1 E^ - (C x /J) j E^ _1 ) + Dq j (33) 

where 6^ and are central-differenced, and the fourth-order smoothing term D is defined by 

Dq j = « E A J (J- 1 ) J |(v n 4 n ) 2 (ji) j + (w e * e ) 2 (Jq) J l (34) 

Here, c^-, the explicit smoothing coefficient must be less than 1/16 for stability. Additionally, 
o=0 for Euler implicit and o = 1/3 for a three-point backward difference. 

An approximately factored form of Eq. (33), which retains the same order of accuracy in t, can be 
obtained if 

[A | + (1 - ajAt^BjmA 3 )- 1 ^ + (1 - oJal^C 3 - ^ « t fi j )|Aq j = LHS(33) + 0 (al) 3 (35) 

(Note: A'* can degrade the factorization error if u is sufficiently small.) Upon substituting the 

left-hand side of Eq. (33) by the left-hand side of Eq. (35), one obtains the factored algorithm. The 
algorithm is solved by the sequence of implicit inversions. 
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[A^ + (1 - a)it(« n B j )lftq« * RHS(33) (36a) 

[A J + (1 - a)a«(^C j - ^ 4 ( fl J )]aq J « A 3 aq* (36b) 

Adding Implicit smoothing to Eq. (36a) 

{*2(1 - til-J]ii(^) n J j+1 D ♦ (1 - a)ac(« n S J )H* = RHS(33) + * I IJ]i 1 («) I 1 J J+1 - •lj 1 ( WA ) n *lj iQj (37) 

where if ej / 0, then e Is usually 2 or 3 times tr and then c E in the right-hand side of Eq. (33) 

Is not restricted. Otherwise If ej = 0, then £ E < 1/16 by stability theory. 

5. GENERATION OF COMPUTATIONAL GRIDS 

The grid-generation process must achieve several ends that Include: (1) developing accurate surface 

representations; (2) distributing points on the body surface; and (3) generating a clustered, well- 
ordered, smoothly varying. Interior mesh. The configurations of current Interest have only moderate vari- 
ations In the axial direction, which simplifies problems associated with grid generation. In general, 
surface-conforming grids are required to simplify application of boundary conditions and to minimize the 
errors Inherent In the thin-layer approximation. Such errors will be kept to a minimum by having the 
n = constant lines coincide with the important velocity and temperature gradients occurring within the 
flow. It is especially Important to make sure the n = constant lines do not Intercept the body at highly 
skewed angles. 

The grid-generation process Is started by choosing c along the body axis. In fact, t Is permitted 
to align with any direction roughly parallel to the body axis. For many problems, a preferred direction 
would be an axis parallel to the lower surface. The t = constant plane, In which n and c vary, is 
chosen as a plane normal to the body axis. Cartesian Inertial coordinates are selected as shown In Fig. 1 
such that x and (, are coincident. This coordinate system simplifies the governing equations since ty 
and = 0. Of more Importance, surface representation and grid generation are reduced to a sequence of 
two-dimensional problems. Continuity of the surface must be maintained from one cross section to the 
next, but this Is usually straightforward. 

The next step Is to represent the surface quantitatively and distribute a surface grid along It. 
Accurate surface representations take advantage of the two-dimensionality of the equations. One simply 
represents each cross section Independently, l.e., with digitized y,z-data points. Intermediate data 
points can be obtained by simple Interpolation. Distributing grid points on the surface Is also greatly 
simplified for the chosen coordinate system because each cross-sectional boundary curve represents a 
simple one-dimensional problem. Previously, It was necessary to distribute points along such curves based 
on prior experiences as to how the solution would develop. 

Finally, a flow-field grid Is generated. This can be accomplished with methods employing either 
algebraic equations or partial-differential equations. The algebraic methods will always require less 
computer-run time and are adequate for a large variety of shapes. The methods employing partial- 
differential equations are preferred for complex shapes to (1) make sure that the n = constant lines are 
(nearly) orthogonal to the surface; (2) allow clustering of grid points into Isolated regions of high 
curvature; and (3) simultaneously satisfy constraints (Implicit In the numerical method) requiring smooth- 
ness and slow variation of the Oacoblans and metrics. 

It should be made clear that a truly orthogonal grid at the body can be (and is) generated only when 
the body shape does not change with axial distance, because the flow-field grid Is confined to a plane 
which Is normal to the axis. This is a direct consequence of the two-dimensional grlddlng methods used 
throughout the prediction procedure. In this context, caution should be exercised In applying the predic- 
tion procedure to bodies having rapid variation in the axial direction. 


6. RESULTS 

A number of different results are presented which show the diversity of the PNS code. The flow 
regimes vary from a Mach number of 2 up to 25. Both laminar and turbulent flows have been considered, as 
well as varying angles of attack. Configurations vary from simple cone-type bodies to lifting winged 
bodies such as the Space shuttle or the generic supersonic cruise fighter. 

The present code was used to calculate the laminar flow field around an ogive-cylinder configura- 
tion. This geometry Is shown in Fig. 2. The cylinder portion begins at an axial location of 
X/L = 0.677 and has a diameter, D. The body origin Is at x = 0. The following are the freestream 
conditions: M,. = 25, Re = 17,000/In., a = 0, the density Is 0.00142 kg/m 3 , and the temperature Is 
268 K. The flow conditions are well within the hypersonic range of Mach numbers and an accurate account 
of real gas effects Is necessary. 
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The grid used has 120 points in the radial direction, 19 points in the meridional direction, and 
518 planes in the axial direction. The marching step size was exponentially stretched from the initial 
value of 8.33 * 10'° to 3.33 * 10 , and then held constant after an axial location of X/L = 0.1667. 

In Fig. 2 the shock position is shown along with the body. The ratio of the shock layer thickness to 
the body radius is only 16* near the nose and increases downstream. This is caused by the expansion over 
the convex curvature of the ogive geometry. At X/L = 0.1667, this ratio is 35* and on the cylinder por- 
tion of the body at X/L = 0.95, this ratio is 115*. 

The second case to be considered is the laminar hypersonic flow over a 30* blunt cone with a 5.6° 
cone half-angle. The free stream Mach number is 14.2; the Reynolds number is 0.62 x 10°/ft; the ratio of 
wall temperature to the stagnation temperature is 0.29; and the angle of attack is 18°. 

Figure 3 shows a comparison of the local separation line as obtained from the experiments (Ref. 17) 
and calculations. The experimental separation line location was determined using the oil-flow tech- 
nique. Two attachment lines (defined as the lines from which the skin-friction lines diverge) coincide 
with the leeward and the windward meridians. In the experiment there was no evidence of any singular 
points at which the axial shear stress is 2 ero. It should be noted that the matrix A is singular for 
zero axial velocity, and the calculations would diverge if there were a reversal in axial velocity. It 
was found that the total shear stress reaches a local minimum along the separation line. Figure 3 also 
shows that the shear stress lines are directed toward the lee in the initial portion of the flow. The 
region of separation starts when the lateral component of the skin friction changes sign in the lee 
region. Both the calculation and the experiment show that the envelope of the skin friction lines exhib- 
its considerable curvature followed by a gradual reduction until the lines finally follow a conical gen- 
erator when the flow becomes nearly conical farther downstream. Also, the skin friction lines, which 
converge around the separation line, emanate from the same region, and the lines near the leeside traveled 
through the open portion of the separated region. 

A systematic study was undertaken to verify the capability of the present computer code in calculat- 
ing the flow field surrounding the flaps on the aft end of the MRVs. 

The configuration studied is a 14“/7° biconic vehicle as shown in Fig. 4, with nose blunting 
(R n = 0.5) and slices on the windward and leeward sides. The control surface protrudes from the cut on 
the aft end of the windward side. Results are presented for M m = 10, a = 10°, T wall = 560° R, and turbu- 
lent flow with Re « 304,800/m. The control surface has a maximum deflection of 15°. 

A typical computational mesh consisting of 45 meridional points by 30 radially stretched points is 
shown In Fig. 5. As can be expected from using straight line rays, a less-than-satisfactory grid struc- 
ture in the region of the flap-flat intersection exists. Using an elliptic grid generator would alleviate 
this problem. Also, in these figures the outermost circumferential line is the fitted bow shock, which is 
tracked by the PNS code. 

A necessary comparison, the pitot-pressure variation through the shock layer at a specified 
x-station on the vehicle, is presented for experimental, inviscid, and laminar flow results. Figure 6 
presents results on the lee side which is just after the biconic juncture. Good agreement between the 
calculated laminar results and the experimental data are apparent. Once outside the viscous layer, the 
Inviscid results also agree quite well. It is seen that the location of the bow shock (as indicated by 
the abrupt drop in the pitot pressure) is predicted quite accurately by both the inviscid and viscous 
computations. These results show the capability of the present PNS code to predict correctly not only the 
surface variables, but also the flow variables Inside the viscous and inviscid portions of the shock 
layer. 

The windward (axial) centerline surface-pressure distributions from just upstream of the flap are 
presented in Fig. 7 for various flap angles. The flap angles are measured from the horizontal and range 
from -7°, which corresponds to the cut (no-flap deflection) to a maximum deflection of 15°, which corre- 
sponds to the maximum angle for no-axial flow separation attainable with the PNS code for the above 
turbulent-flow conditions. Note the rapid compression at the corner of the flap which signifies the 
location of the flap shock, followed by a lesser compression as the flow continues up the flap. These 
results are quantitatively correct when compared to the nonseparated experimental results of Settles 
et al. (Ref. 18) and also are consistent with the results presented in Ref. 19 which were for a slightly 
different Mach number and configuration. 

Figure 8 presents a comparison between the experimental and calculated heat-transfer rates on the 
4 * = 90° and the « = 180° (lee side) meridians. The slight discrepancies are most likely because of a 
problem with predicting flow transition between laminar and turbulent flow. 

The PNS code has been used to determine the viscous supersonic flow over complex configurations. One 
such configuration geometry consisted of a thick-finned projectile (Fig. 9). The free-stream conditions 
were M^ = 4, a turbulent boundary layer with Re/L = 10°/1n., T^ = 100°R, T w _ i 1 = 540°R, and a « 0° 
and 2°. Forty-five grid points were stretched In the radial direction, and 117 grid points were distrib- 
uted in the circumferential direction. The starting solution at x= 2.0 In. was generated with the PNS 
code using the conical step-back procedure described earlier. The conical starting solution was then 
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inarched over the cone-cylinder portion of the projectile to x = 13.55 in. (which is just in front of the 
fins). 

The grids used to calculate the flow over the finned portion of the projectile were generated using 
an elliptic grid generator similar to the one developed in Ref. 20. This procedure furnishes both spacing 

and angular control of the grid lines as they approach the boundaries of the domain such as the body and 

the shock. Figures 10a and 10b show representative grids at various streamwise x locations. The 
n » constant grid lines Intersect the body orthogonally. This is Important for the accurate Implementa- 
tion of the viscous surface boundary conditions. The solution at x = 13.55 in. was then marched to the 
end of the projectile, using the type of grids shown in Figs. 10a and 10b. 

Figure 11 shows the density contours obtained for a six fin case. The protrusion of the fins into 
the supersonic flow field causes fin shocks, shown by the clustered contours above the fins. The clus- 
tered contour levels close to the body surface, however, are caused by the large heat-transfer rates that 
occur because of the constant wall -temperature boundary condition that has been Imposed. 

In Table 1, a comparison between numerical and experimental force coefficients (Sturek, W. B. , 
private communication, 1982) is presented for the six-finned projectile. 

The numerical and experimental values of C N a and x c -p differ by less than 4%. In comparing the 
values of C A , It should be kept In mind that the base drag has not been added to the numerically obtained 
value of C A . It should also be noted that there are differences between the projectile geometry used In 
the PNS calculation and in the experiment. The experimental projectile has a blunt nose and grooves on 

the cylindrical portion of the forebody In order that It could fit Into a sabot. . 

Another complex configuration Is the Space Shuttle orbiter. Numerical results have been obtained for 
the following wind-tunnel conditions: M„ = 7.9, a = 25°, T wal1 = 540°R, Re » 60,728/in. turbulent flow. 
For this calculation, the Shuttle surface coordinates were obtained from Rockwell-International Corpora- 
tion. The current geometry consists of the complete Shuttle: the canopy, 0MS pods, and the vertical 

stabilizer are Included. 

The three-dimensional blunt-body code originally developed by Kutler et al. (Ref. 15) was used to 
obtain the blunt-nose solution which creates the necessary starting planes for the PNS code at 
X/L * 0.0522. This solution was then marched downstream using the elliptic grid generator to construct 
the grid between the body and the fitted outermost shock wave. The grid consisted of either 61 or 121 
points in the meridional direction, and 45 geometrically stretched radial points. An example of the grid 
at an X/L = 0.66 Is shown In Fig. 12. The outermost grid line Is the bow wave, which is fitted using an 
Implicit technique. 

The computer-generated particle paths of Fig. 13 exhibit the main features of the flow field which 
surrounds the Shuttle. Two distinct phenomena are the vortices on the leeslde of the body and the vortlce 
owing to the strake wing. 

Finally, a generic supercruise fighter configuration is considered. Numerical results for supersonic 
cruise at M m = 2.169 are presented. The wind-tunnel conditions considered are such that the Reynolds 
number is turbulent. Adiabatic wall conditions are assumed at the body's surface. Angle of attack of 4° 
is considered. The absence of yaw results in a pitch plane of symmetry which reduces the computational 
space and, hence, the amount of central-processing-unit (CPU) time and storage requirements. The current 
geometry, consisting of a canopy, body, nacelle and wing. Is In a form that Is readily usable by the PNS 
code. The point distribution on each cross section was determined Independently and was then used to 
produce a geometry file which was read into the PNS code. Intermediate cross sections were determined by 
simple Interpolation in the axial direction. All of the figures are scaled by either the model length (L) 
(Fig. 14), or the local span width (b). The locations of some of the referenced cross sections are shown 
In Fig. 14. 

A typical grid, obtained by using the elliptic mesh generator, is presented in Fig. 15. The outer- 
most boundary Is the location of the fitted bow shock. This particular grid has 61 meridional points and 
45 radial points. 

Figures 16 and 17 show mesh studies for two different X/L stations. Each study used 45 radial 
points, but varied the number of meridional points: 61, 91, and 121 for o = 4°. The results compared 
well with experimental results. The coarse details of the flow field were reasonably resolved by all 
systems; however, for the X/L = 0.8 station it appears that the 121 * 45 grid system produces the more 
accurate overall results. . . 

Recently, a numerical result for = 2.2 flow past a wlng-fuselage-nacelle has been completed. In 
Fig. 18, Cp distributions for some representative cross sections are presented. The solid line corre- 
sponds to a wing-fuselage result, the open symbols correspond to a wing-fuse lage-nacelle result, and a 
solid symbol corresponds to experimental results. Additionally, each open symbol Is the surface location 
of the grid points used In the calculation. As expected, the upper surface Cp distribution Is not 
affected by the nacelle on the under surface of the wing. The Cp distribution In the vicinity of the 
nacelle Is qualitatively correct and, for the sparse experimental data, shows a good comparison. 
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7. CONCLUOING REMARKS 


A PNS code has been presented which was used to solve the high-speed flow about a wide variety of 
complex configurations. The code which has evolved over the last eight years is a viable tool to investi- 
gate the aerodynamics and the fluid physics of supersonic flow past complex configurations. Most 
recently, an additional capability of the PNS code has been developed. This consists of using nonbody 
axis normal planes to march downstream over complex configurations. At the trailing edges of lifting 
surfaces, this procedure might be the best since one could orient the marching planes as to march off the 
surfaces parallel to the trailing edges. 
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Table 1 

Comparison of experimental and numerical force coefficients. 


Force 

Lift curve slope, 

Location of center of pressure. 

Axial force coefficient, 

coefficient 

C N (per radian) 

a 

X C.P. (calibers) 

C A 

Numerical 

10.395 

9.361 

0.304 

Experimental 

10.750 

9.745 

0.297 
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Fig. 6. Comparison between computations and 
experiment of the pitot-pressure variation through 
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Fig. 13. Computational particle paths up to 
X/L = 0.66. 
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Fig. 18. Cp variations at different cross sections for wing-fuselage-nacelle at M„ = 2.2 and a = 10°. 
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